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SOUND  FIELD  COMPUTATION  BY  MEANS 
OF  DESK  TOP  COMPUTERS 


By 

Audun  Skretting 


ABSTRACT 


Several  formulae  concerning  sound  in  sea  wafer  are  presented  in 
forms  that  are  convenient  to  program  on  a  desk-top  computer  (DTC). 

They  are: 

a)  Formula  for  calculating  sound  speed  in  sea  water. 

b)  Ray  tracing  formulae  for  calculating  the  sound  ray 
path  in  a  medium  with  the  sound  speed  varying  only 
as  a  function  of  depth. 

c)  Intensity  formulae  for  calculating  the  intensity 
along  a  ray  already  plotted. 

d)  Travel-time  formulae  for  calculating  the  travel-time 
for  a  ray  already  plotted. 

e)  Pressure-depth  formula. 

Special  applications  to  one  existing  DTC  (the  Olivetti  101)  are 
given,  with  the  programs  given  in  tables.  Two  examples  to  illustrate 
the  program* s  usefulness  are  given  at  the  end. 
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INTRODUCTION 


The  computation  of  sound  fields  in  underwater  acoustics  started 
mainly  during  World  War  II  i_Ref«  l]» 

From  then  to  the  late  1950 1 s,  when  ray-tracing  programs  for 
digital  computers  became  more  and  more  common,  many  different 
approaches  were  made  to  the  problem  I^Refs.  2  and  3]*  The 
numerical  solution  is  not  difficult,  but  to  give  a  picture  of 
the  sound  field  requires  many  calculations  to  be  made  with 
far  greater  accuracy  than  is  possible  with  a  slide  rule* 

The  solutions  were  usually  obtained  with  mechanical  calculators, 
analogue  computers,  or  numerical  calculations  by  means  of 
precomputed  tables* 

As  access  to  big  computers  became  common,  these  took  over  the 
task  of  ray-tracing. 

During  the  last  years,  different  types  of  desk-top  computers 
(DTC)  have  been  put  on  the  market.  They  are  digital 
"mini-computers11  ,  and  despite  their  low  price  they  have 
memory  for  both  instructions  and  figures.  This  paper 
demonstrates  how  these  can  be  used  to  solve  ray-tracing 
problems . 

The  programs  presented  here  for  desk-top  computers  are  not 
intended  to  compete  with  ray-tracing  on  large  computerso  But 
when  a  large  computer  is  not  available  (on  a  research  vessel, 
say,  or  in  a  small  factory)  or  when  it  is  overloaded,  they 
are  quite  helpful  for  small  ray-tracing  problemso  Also  a 
short  computation  of  the  sound  field  by  a  DTC  can  be  helpful 
when  specifying  the  data  for  a  large  computer's  ray-tracing 
program. 
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The  report  is  divided  into  three  parts. 


Chapter  1  presents  formulae  for  computing  the  sound  field, 
with  the  simplifications  necessary  for  a  DTC. 

Chapter  2  contains  the  application  of  the  formulae  for  an 
existing  DTC  (the  Olivetti  101). 

Chapter  3  presents  two  examples  of  ray-tracing,  illustrating 
the  use  of  the  different  programs. 

The  appendix  contains  instructions  for  use  of  the  given 
programs  on  an  Olivetti  101  DTC.  It  also  contains  a  short 
program  with  full  description  of  the  program  steps  so  as  to 
provide  an  insight  into  the  programming  methods. 
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1 


THE  FORMULAE 


1  •  1  Sound  Speed  Formula 

Wilson* s  second  formula  is  at  present  widely  used  to  calculate 
the  speed  of  sound  in  sea  water0  Its  complexity,  however, 
necessitates  the  availability  of  electronic  computers,  for 
which  it  was  primarily  establishedo 

In  Ref «  4 >  C®C.  Leroy  proposed  two  formulae:  one  that 

approaches  Wilson's  second  formula,  and  one  that  approaches  the 
data  Wilson  used  when  computing  his  second  formulae  Although 
the  formulae  proposed  by  Leroy  are  simple,  they  are  fairly 
accurate  and  it  is  demonstrated  that  the  latter  formula  fits 
Wilson's  data  better  than  Wilson's  equation0  The  two 
formulae  are  rather  similar,  and  here  the  second  is  dealt  with. 

The  formula  is  presented  in  Table  1.  It  has  been  divided  into 

a)  A  minimum  number  of  terms,  VQ,  forming  the 
’feimplif ied"  formula,  which  is  sufficient 

with  a  slightly  reduced  accuracy  for  depths  of 
less  than  1000  m  at  any  latitude  provided  that 
t  <  2  5°  C  and  3(%)<  s  <  42<&,  . 

b)  A  depth,  temperature,  and  latitude  dependent 
expression,  V_^,  which  is  necessary  to  operate 
deeper  than  1000  m.  It  extends  the  validity  of 
the  formula  to  depths  of  7000  m  and  improves 
the  accuracy  to  ±  0.1  m/s  (for  t^23°C)o 

c)  A  temperature  correction,  V^,  for  t  ^  2  5°  Co  This 
correction  still  leaves  the  results  at  lower 
temperature  within  0o 1  m/s  of  Wilson's  formula 
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and  permits  the  calculation  of  sound  speed 
for  temperatures  up  to  34° Co 

d)  A  correction  term,  ,  for  very  great 

depths,  necessary  only  if  d  >  7000  md 

e)  A  corrective  term,  V^,  for  low  salinity 
to  be  applied  only  when  s  <  30&,  • 


It  follows  that  the  way  to  make  programs  for  sound  speed  when 
using  Leroy1 s  formula  is  to  make  one  program  for  the  basic 
formula,  and  two  small  programs  for  the  deep-ocean  term  and 
the  low-salinity  termc 


1  *  2  Ray  Tracing  Formulae 

Figure  1  shows  a  sound  ray  in  a  layer  of  thickness  Az  with 
constant  sound  speed  gradient* 


5 


TABLE  1 


FORMULA  FOR  THE  CALCULATION  OF  SOUND  SPEED  IN  SEA  WATER 


COMPLETE  V  =  V„  + V  +  V,  +  V  +  Vj 

-  0  a  b  c  d 

BASIC  V  =  Vn  +  V  +  V, 

-  0  a  b 

SIMPLIFIED  V  =  V0 

in  which 

v0  =  1493  +  3 ( T  -  10)  -6  X  10"3(T  -  10)2  -4  X  10"8  (T  -  l8)2 
+  1.2  ( S  -  3  5)  -10"2(T  -  18)(S  -  35)  +z/6i 

V  =+10-1C2  +2  x  10_4C2  (T  -  18)2  +10-1C$/90 

3 

Vb  =  2.6  x  10"4t(t  -  5) ( t  -  25) 

vc  =  -10_3C2(C  -  4) ( C  -  8) 

Vd  =  +l'5  *  10~3(s  -  3  5)2(  l  -  c)  +3  x  l0_6t2(t  -  30) (S  -  35) 


Where 


V  is  the  sound  speed  in  m/s 
T  is  the  temperature  in  °C 
S  is  the  salinity  in 

z  is  the  depth  in  m,  and  Q  =  z/l000?  the  depth  in  km 
$  is  the  latitude  in  degrees 


Note :  VQ  can  also  be  written: 

VQ  =  1449.44  +  4 • 56T  -  0.046  T2 

+  1 . 2( S  -  3  5)  -10-2(T  -  18) (S  -  35)  +z/61 
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From  Snell Ts  law  we  have 


n 


m 


cos  0 


n 


cos  0™ 
m 


k  , 


where  k  is  a  constant  for  the  ray  in  any  layer. 

We  introduce 


[Eq.  l] 


Ac  =  c  -  c 
m  n 


Ref.  5  gives  the  horizontal  distance  in  one  layer J 


A  v  A  z  •  k  /  .  X 

AX  =  •  (sin  en-sm  0m> . 


As  we  have  from  Eq®  1  : 

c 

n 

cos  9  —  -j—  , 

un  k  9 

we  get 


/  /  / cn\ 2 
sin  en  =  y  i  -  coS2  en  =  j  i  -  j 


[Eq.  2] 


[Eq.  3] 


[Eq.  4] 


[Eq.  5] 


By  introducing  Eq.  5  into  Eq.  3  we  have  an  expression 
for  the  horizontal  distance  for  which  it  is  possible  to  make  a 
program.  But  if  we  examine  Eq.  3>  we  notice  that  it  contains 
two  terms: 


A  z  • 


Ac 


k 


[Eq.  0] 
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where  r  is  the  radius  of  the  circle  of  which  the  sound  ray 
forms  an  arc  in  the  layer  ,  and  the  term 


[Eq.  7] 


sin  0  -  sin  g  , 

un  um  7 


If  we  have  a  layer  of  nearly  isovelocity  water,  Eq.  6  will 
have  a  very  large  value  and  Eq.  7  a  very  small  value  and  the 
product  will  be  rather  unreliable  because  the  computer  must 
work  with  a  restricted  number  of  decimals.  To  alter  the 
relative  size  of  the  two  terms,  and  also  to  have  a  formula 
that  is  easy  to  evaluate,  we  use  an  approximation  for 


From  Eq.  5  we  have 


As  k  +  c  ~  2k,  we  get 


sin  Qn~  k 


1 


k  -  c 


n 


/ 


k  -  c 


n 


[Eq.  8] 


Equations  8  and  3  give 


By  comparing  Eq.  3  and  Eq.  9  we  see  that  the 


r 


[Eq.  10] 
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A 


of  this  formula  is  reduced  about  27  times  compared  to  r  [Eq*  6]. 
This  means  that  the  term 


[Eq.  1 1] 


is  also  increased  27  times  compared  to  Eq.  7<> 

We  still  have  the  problem  that  if  Ac  goes  to  zero  for  one 
layer,  Eq*  10  goes  to  infinity  and  Eq.  11  goes  to  zero. 

But  this  can  easily  be  arranged  in  the  following  way. 

The  sound  speeds  (in  m/s)  put  into  the  computer  need  be  given 
to  two  decimal  places  only,  because  as  the  sound  speed  formula 
and  the  sound  speed  meter  have  a  standard  error  of  about 
0.3  (m/s)  there  is  no  justification  for  taking  more  decimal 
places*  The  value  of  Ac  must  go  in  steps  of  O.Olo  If  Ac  =  0, 
it  can  be  made  equal  to  0.01  in  the  computer,  and  the  error 
introduced  is  surely  less  than  the  error  made  by  uncertainty 
in  sound  speed* 

In  Eq.  9,  Ad  must  always  be  positive  and  Ac  must  always 
be  equal  to  c  -  c  to  obtain  a  positive  value  for  AX* 

Equation  9  is  rather  simple  to  compute.  Some  precautions 

must  be  taken  however.  If  we  give  the  computer  a  sound  speed 

bigger  than  k,  we  may  still  get  a  result  because  some  computers 

take  the  square  root  of  the  absolute  value,  thereby  ignoring 

the  fact  that  in  this  case  we  have  a  negative  number  and  thus 

cn 

have  no  solution.  But  as  k  =  -  is  constant  for 

cos  0 
u  n 

the  ray,  the  ray  is  not  able  to  penetrate  into  an  area  with 
sound  speed  greater  than  the  k  value.  In  fact  the  ray  will 
turn  when  it  reaches  a  sound  speed  equal  to  its  k  value, 
because 


k 


c 


n 

cos  0 
un 


c 

n 


cos  0 


c  for  0  =  0°  * 

n  vn 
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k 


That  means  that  in  the  program  we  must  test  whether  the  given  cn 
is  bigger  than  the  k  value  and?  if  so,  compute  the  turning 
point  and  modify  the  input  datao 

The  computation  of  the  turning  point  is  simple  because,  as 
for  the  rest  of  the  computation,  we  suppose  linear  gradient. 

The  depth  of  the  turning  point  (where  c  =  k)  is 


Zturn 


z 

n 


[Eq.  12] 


c 

n 

In  the  program  we  need  the  value  of  k  =  Fos — 5”  >  n  being 

°n 

any  point  on  the  ray.  Normally  k  is  either  determined 
by  the  initial  angle  and  sound  speed,  or  by  a  turning  point 
(where  k  =  c^ )  •  To  compute  k  from  the  initial  data  a 
separate  program  must  be  made  where  cQ  and  0^  are  given 
as  input So 


Cosine  0  can  be  computed  by  means  of  the  series 


cos  0  = 


1 


X2  ,  X4 
~  2  41  ’ 


[Eq.  13] 


where  X  is  the  angle  in  radians.  By  using  only  three  terms, 
one  makes  an  error  at  30°  of  about  3  x  10  ,  which 

corresponds  to  less  than  0.1°.  As  the  ray-tracing  program 
should  not  be  used  with  angles  bigger  than  about  30°  because 
of  the  simplification  in  the  formula,  this  accuracy  for  the 
cosine  function  should  be  sufficient. 
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1 o  3  Intensity  Formulae 


Reference  .5  gives  the  formulae 


T~ 

p 


sin  ep 

c°s  0O 


dX 


[Eq.  14] 


which  is  valid  regardless  of  the  path  followed  by  the  ray 
between  the  starting  point  and  the  point  of  interest  p. 

To  compute  this  equation,  we  need  an  expression  for  - -  • 

.  .  .  .  0 

Assuming,  as  before,  that  we  can  divide  the  ocean  into 

horizontal  layers  of  constant  sound- speed  gradients,  Ref o  5 

gives  the  formula: 


n 


w 


dX 

d0c 


n 


m 


=  r  tan 

n 


0o( 


l 


sin  9 


m 


Sln  9n 


[Eq.  15] 


Here,  r  is  the  curvature  radius  in  layer  n,  and  a  is 
n  dX  n 

the  amount  by  which  3-5 —  is  changed  by  passing  from  layer  m 


to  layer  n.  Therefore: 


dX 

_ LL 

de0 


p 

s 


n=  1 


a 


n 


[Eq.  16] 


By  using  the  same  approximation  for  sin  0n  in  Eq.  15 
as  before  [see  Eq.  8]  we  get: 


a  =  Az 

n  Ac 


V  2k  -  tan  Q( 


±  1  ±  1 


Vk~ 


c 

m 


[Eq.  17] 


The  ray-tracing  program  must  print  out  the  value  [Eq.  10]: 


r 


t 

n 


72  k 


1 1 


Az 

Ac 


and  it  will  be  used  as  input  data  here 


Combining  Eqso  10,  14*  16  and  17  gives: 


sin 


tan  0O 

.  I  - - 

n  cos  0 


±J—  _  .  ±  1 


[Eq0  18] 


In  special  cases  there  are  some  problems  in  evaluating  the 
formula;  these  are  all  discussed  in  Ref.  5» 

The  first  problem  is  when  0Q  =  0o  Then  Eq<>  6  cannot  be 
used,  and  with  the  program  made  [see  Ch.  2]  it  is  not  possible 
to  calculate  the  intensity  for  a  ray  with  an  initial  angle  of 
zero. 

The  next  problem  is  when  the  ray  turns  at  point  n  so  that  0n  =  ^° 
If  we  look  at  Eqe  15*  we  see  that  it  is  not  possible  to  use 
this  point.  But  we  have  r  =  r  ,  -  and  sin  fl  , „  =  -sin  Q  . 

Thus  we  can  drop  point  n  and  go  to  n+1.  Then  sin  Q  will 
change  sign,  and  in  the  computation  of  Eqo  18  we  must  change 
the  sign  of  the  approximated  sin  9  term.  This  can  be  done 
by  storing  a  ±1  to  use  in  the  denominator,  and  changing  its 
sign  every  time  a  turning  point  is  passed.  By  this  simple 
approach,  it  is  not  possible  to  get  the  intensity  in  the 
turning  point. 

1 o  4  Travel- time  Program 

When  a  ray  tracing  is  made,  one  often  realizes  that  a  certain 
point  can  be  reached  by  several  sound  rays.  To  find  the 
succession  of  the  different  arrivals,  it  is  necessary  to  compute 
the  travel  time. 
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From  Ref o  6  we  have  the  formulae 


t 


k2 


dz 


[ Eq.  19] 


for  the  travel-time  in  seconds 


X  = 


1 

k 


U 


dz 


[Eq.  20] 


Differentiating  Eq.  20  gives 


dX 


[ Eq.  2 1 ] 


Eliminating  dz  from  Eq.  19  by  means  of  Eq.  21  gives 


t  =  $  +  k 
k 


r  (.L . 

J  \  c2  k 2  / 


dX 


[Eq.  22] 


(k+  c)  (k  -  c) 

k2  c2 


dx«$  +  -2- 

k  k» 


(k  -  c)  dX 


[Eq.  23] 


Here  we  have  made  the  approximation  that  k  4-  c  2k  and 
k2  c2  k4  • 
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The  integral 


(k-c)dX  is  not  easy  to  evaluate,  since  c 


is  a  linear  function  with  depth  in  each  layer,  but  is  a 
complex  function  of  the  horizontal  distance  X.  However, 
by  selecting  the  layer  thickness  so  that  sound  speed 
differences  are  not  too  great  within  each  layer,  we  can 
use  the  approximation 


c(  a)  c  + 
m 


a  . 


[Eq.  24] 


Equations  5  and  6  give 


X 

-2  + 

k 


2  Ax 


n=l 


c  +  c 
m  n 


[Eq.  25] 


Equation  25  contains  the  k  value,  the  sound  speed  of  each 
layer,  and  the  horizontal  distance.  It  will  normally  be 
difficult  to  find  space  in  a  DTC  to  calculate  the  travel-time 
in  the  main  ray-tracing  program.  But  when  the  ray-tracing 
is  done  and  the  horizontal  distance  Xp  is  known,  a  small 
and  easy  program  finds  the  travel-time. 
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APPLICATION  OF  THE  FORMULAE 


2  . 


2.1  Introduction 


There  are  already  a  number  of  different  DTC 1 s  on  the  market  and 
the  number  will  surely  increase  in  the  future.  To  give 
recommendations  for  the  programming  of  different  sizes  of  DTC 
is  outside  the  scope  of  this  report.  This  section  gives  and 
explains  the  flowcharts  that  are  linked  to  the  special  DTC  for 
which  the  programs  are  made  (Olivetti  101),  but  they  are 
general  and  can  serve  as  guides  when  programming  on  other  DTC’s. 
However,  if  the  DTC  being  used  has  smaller  memory  and 
programming  facilities,  the  programming  of  the  presented 
formulae  might  be  difficult  or  impossible. 

If  some  details  of  the  explanation  are  not  clear,  it  might  be 
helpful  to  study  the  Appendix , which  gives  a  short  description 
of  the  actual  DTC,  the  programs  in  table  form,  and  a 
description  of  the  use  of  the  programs.  The  Appendix  also 
contains  an  example  in  which  each  step  of  the  programming 
is  explained. 

The  numbers  indicated  on  the  flowchart  [Figs.  2,  3,  4  and  5] 
correspond  to  instruction  numbers  in  the  corresponding 
tables  [Tables  9*  12,  13>  14  respectively]."*  The  number  of 

constants  that  can  be  put  into  a  register  and  the  number  of 
constants  it  is  necessary  to  generate  in  the  computer  will 
vary  from  computer  to  computer.  The  actual  organization  for 
the  Olivetti  101  is  therefore  only  mentioned  in  the  Appendix. 

The  programs  for  computation  of  the  soundfield  are  all  based 
on  a  sound-speed  distribution  varying  only  as  a  function  of 
depth.  The  first  step  is  therefore  to  approximate  the  sound 


See  programs  between  coloured  pages  at  end  of  main  text. 
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speed  profile  with  straight  segments,  thus  dividing  the  medium 
into  layers  of  constant  sound- speed  gradient «  If  data  are 
needed  at  a  depth  that  does  not  coincide  with  a  depth 
where  the  gradient  is  changing,  it  is  necessary  to  compute 
the  sound  speed  at  this  depth  This  will  often  be  the  case 
for  the  source  depth. 

The  programs  work  in  series s  one  first  uses  the  sound-speed 
program  to  get  the  sound  speed  profile  and  then  uses  the 
ray-tracing  program.  These  provide  input  for  the  intensity 
program  or  the  travel-time  program. 


2.2  Sound  speed  program  [Fig.  2 

The  formulae  presented  in  Sect.  1.1  can  easily  be  matched  to 
different  sizes  of  DTC.  Here  the  formulae  are  programmed  in 
two  programs 

a)  One  program  for  the  basic  formula  V  =  Vo  4*  4-  V^- 

b)  One  program  contains  the  great-depth  term  Vc 

and  the  low  salinity  term  Vd .  They  are  programmed 
in  series,  so  one  calls  the  program  for  Vd  or 
Vc  by  using  different  labels  on  the  keyboard. 

The  program  for  the  basic  formula  will  cover  most  ocean  cases. 
When  the  depth  is  greater  than  7000  m  or  the  salinity  is  less 
than  30 %r  ,  one  has  to  use  the  other  program  to  find  the 
correction  term. 

The  term  Vo  is  rewritten  to 


Vo  =  1401.14  +  49.1  t/10  -  4-6  (t/10)2  4-  13.8  s/10  -(t/l0)(s/l0)  +  z/61  . 


Only  the  flowchart  for  the  main  program  is  shown  in  Fig,  2,  the 
other  program  is  so  short  and  simple  that  a  flowchart  would  be 
superfluous , 
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2,3  Ray  tracing  program  [Fig*  3] 


Because  of  space  problems,  one  short  program  (not  illustrated) 
is  used  to  compute  the  k  value: 


cos  e0 


For  the  DTC  with  cos  0  or  sin  0  as  a  built-in  subroutine, 
this  is  of  course  very  simple.  Here  cos  0  is  computed  by 
using  Eq.  13.  The  k  value  is  then  used  as  an  input  parameter 
to  the  ray  tracing  program. 

For  the  ray  tracing  program  Eqs.  9  and  12  are  used.  It  must 
also  be  remembered  that  the  program  must  print  out  tan  0C  and  r! 
[Eq. 10]  for  the  intensity  program  [Sect.  2.4]*  As  tan  0  had  to 
be  computed,  the  program  was  arranged  so  that  this  computation 
and  printing  is  done  for  all  0  . 

The  start  of  a  program  is  reached  by  calling  a  label.  The 
first  part  of  the  program  receives  the  initial  data,  the 
k  value,  cQ  (sound  speed  at  source),  and  z0  (depth  of 
source),  and  arranges  them  in  registers.  Take  note  of  Instr.  8 
which  cleans  Register  D  after  new  initial  data;  the  sum  of  AX 
(horizontal  range  in  each  layer)  is  stored  in  this  register. 

At  the  end  of  the  initial  data  section,  the  computer  is  told 
to  jump  to  the  last  part  of  the  program  for  computing  and 
printing  tan  @c. 

The  main  loop  of  the  program  first  receives  the  sound  speed 
of  the  actual  layer  as  input  data.  It  stores  the  value,  and 
computes  and  stores  the  Ac.  In  Instr.  18  the  program  checks 
if  Ac  =  0,  and  if  so,  the  sound  speed  c  and  the  Ac  is 
modified  by  0.01  m/s. 

At  Instr.  28,  the  depth  of  the  layer  (z  )  is  given  as  input, 
stored,  and  the  Az^  computed. 
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FIG.  3  FLOWCHART  -  RAY  TRACING  PROGRAM 
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Instr.  35  compares  the  sound  speed  with  the  k  value.  If  the 
sound  speed  is  smaller  or  equal  to  the  k  value,  the  computer 
jumps  to  Instr „  55*  If  the  sound  speed  is  greater  than  the 

k  value,  then  we  have  asked  for  a  sound  speed  and  corresponding 
depth  that  the  ray  does  not  reach,  because  it  turns  beforehand. 
The  computer  then  computes  the  turning  point  [Eq.  12!,  prints 
out  the  value,  and  corrects  the  values  for  c  ,  Ac  ,  z 

?  n  ;  ti  J  n 


At  Instr.  56  the  evaluation  of  the  formula  starts.  Instruction  74 
prints  out  the  value  of  r*  [see  Eq.  10],  which  is  needed  as 
input  for  the  intensity  program.  As  we  have  the  same  horizontal 
range  whether  Az  is  positive  or  negative,  we  have  to  take 
the  absolute  value  of  the  result.  The  program  sums  all  the 
results  in  Register  D  and  prints  out  the  sum. 

Instructions  79  be  90  compute  tan  and  print  out  the  value. 

This  figure  is  not  necessary,  only  tan  0O  being  used  for 
the  intensity  program.  By  putting  these  instructions  between 
Instrs.  8  and  10,  the  computer  program  will  be  a  bit  faster. 

But  very  often  it  will  be  convenient  to  know  the  angle  in 
the  various  layers. 

2 . 4  Intensity  program  [Fig.  4 1 

Equation  18  is  used.  The  evaluation  of  the  formula  is  made  under 
the  assumption  that  one  is  only  interested  in  the  intensity  at 
certain  points  on  the  ray.  Only  the  term 


tan  0O 


Y 


Y 


r  i  - 

n  cos  0O 


[Eq.  26] 


(where  y  =  +1  or  -l) 

is  evaluated  and  summed  for  each  set  of  input  data 
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When  the  intensity  at  a  point  is  wanted,  one  calls  the  last  part 
of  the  program  by  pressing  a  label.  Horizontal  distance  X 
is  given  as  an  input  parameter,  and  the  result  is  printed  out. 

For  the  final  result,  the  product  t;  sin  0  is  needed.  It  is  easy 

Z  P 

to  find,  by  means  of  the  same  approximation  for  sin  0  as  used 
before,  that 


k 

2  sin 


[Eq.  27] 


Since  I./l  can  often  be 
o/  p 

not  to  get  overflow  in  the 
difficult  to  read,  e.g.  a 
printed  out  as  31620000. 


printed  out  as  J  T^JT~  = 
pressure  1  metre  from  the 
point  p. 


a  large  number,  care  must  be  taken 
registers.  Large  numbers  are  also 
not  unusual  loss  of  75  dB  will  be 
In  the  program  the  result  is  therefore 
P0/Pp>  where  P^  is  the  initial 
source  and  P  is  the  pressure  on 


This  amount  is  equal  to  the  horizontal  distance  when  we  have 
spherical  spreading,  and,  for  example,  the  same  75  dB  loss  will 

fl  Q  f  1q 

be  printed  out  as  5624-  The  value  of  - a^so  Printed 

out  to  show  how  the  loss  compares  to  normal  spherical  spA  eading. 

Instructions  1  to  19  take  care  of  the  initial  data,  and  the 
main  program  loop  starts  at  Instr.  20. 

Instructions  24  to  34  check  if  c^  —  c^,  and  if  so  modify  c^. 
This  is  to  work  in  accordance  with  the  ray  tracing  program. 

Instructions  35  to  44  check  if  c^  =  k,  and  if  so  the  sign  of 
the  denominator  ( y)  in  the  approximated  sine  term  is  altered. 
The  computer  then  stops  [instr.  42]  to  wait  for  the  next  input 
data,  and  returns  to  the  start  of  the  main  loop.  Note  that 
nothing  is  done  with  the  last  input  data,  and  Instr.  42  is  only 
to  take  care  of  the  sequence  in  the  input  data. 
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A  computation  of  Eq.  28  starts  at  Instr.  46  and  finishes  at 
Instr.  56.  The  result  is  added  to  and  stored  as  result  in 
the  C  register. 

When  the  intensity  is  wanted,  the  last  part  of  the  program  is 
called  by  means  of  a  label. 

Instructions  65  to  70  compute  Eq.  27  and  the  final  result  is 
printed  out  at  Instr.  74° 

An  initial  value  for  y  =  -hi  has  to  be  put  into  the  D 
register . 

It  is  not  possible  to  calculate  the  intensity  for  a  ray  starting 
horizontally,  or  to  ask  for  the  intensity  at  a  point  where  a  ray 
is  horizontal , 

For  a  splitting  ray,  it  is  possible  to  compute  the  intensity  for 
the  ray  turning  back,  but  not  for  the  ray  bending  away. 

If  one  wants  to  calculate  the  intensity  for  a  reflected  ray, 
one  has  to  shift  the  sign  on  all  the  r^  (given  as  input  data) 
after  each  reflection.  This  means  that  one  gives  the  r^  as 
they  are  printed  out  from  the  ray  tracing  program  up  to  and 
including  the  reflection  point.  Thereafter  one  gives  the  r^ 
with  opposite  signs  to  those  printed  out  until  the  next 
reflection,  and  so  on. 

2.5  Travel -time  program  [Fig.  5] 

Equation  25  is  used.  The  initial  data  are  k  and  cQ.  The  input 
data  are  sound  speed  of  each  layer  the  ray  is  passing, and  horizontal 
distance.  The  program  checks  if  the  sound  speed  is  greater 
than  the  k  value,  and  if  so,  it  substitutes  the  sound  speed 
by  the  k  value. 
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Output  from  the 
anomaly”  At^  = 
convenient  when 


program  is  travel-time  t^  and  "travel-time 
t^  -  Xn/c0.  The  last  figure  is  very 
comparing  travel-times  for  different  rays. 
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ILLUSTRATION  OF  THE  PROGRAMS 


3 . 1  Complete  Ray-tracing  Examples 
The  data  are  taken  from  Ref,  7 • 

Table  2  gives  the  depth,  temperature,  and  salinity  together  with 
the  sound  speeds  computed  by  means  of  the  sound  speed  program. 

As  the  salinity  in  this  case  is  more  than  30%o  and  the  depth 
is  less  than  7000  m,  the  use  of  the  ”basi.c  formula”  is  sufficient. 


TABLE  J? 


MEDITERRANEAN ,  43°23'N 

Date  22  September 

0°47'E 

1956 

Depth 

Temperature 

Salinity 

Computed  Sound  Speea 

0 

21.0 

38.25 

1520.37 

10 

20.4R 

30.24 

1527.13 

20 

19.36 

38.22 

1524.22 

30 

14.69 

38.14 

1511.11 

50 

13.20 

38.20 

1506.10 

75 

13.01 

38.30 

1506.01 

90 

13.09 

38.35 

1506.59 

160 

13.13 

38.40 

1507.93 

210 

13.42 

38.49 

1509.82 

310 

13.42 

38.53 

1511.51 
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Figure  6a  shows  the  sound  speed  profile* 


Suppose  now  that  we  have  a  hull-mounted  sonar  with  a  transducer 
at  10  m  depth*  (As  the  source  depth  then  coincides  with  one 
of  the  points  of  the  sound  speed  profile,  the  computation 
of  the  sound  speed  for  the  source  depth  is  superfluous) * 

The  following  will  explain  how  to  calculate  the  ray  limiting 
the  shadow  zone, and  some  other  rays* 

By  means  of  the  k-  value  program,  one  computes  the  k  value 
for  the  desired  rays,  giving  sound  speed  and  initial  angles 
at  the  source  as  input  parameter.  As  this  shape  of  sound  speed 
profile  gives  a  downward  refraction  of  the  rays,  the  limiting 
ray  will  be  the  one  that  is  horizontal  at  the  sea  surface.  The 
k  value  for  a  ray  is  equal  to  the  sound  speed  at  the  depth 
where  it  is  horizontal,  hence  for  this  ray  the  k  value  is  equal 
to  the  sound  speed  at  the  surface* 

The  rays  from  -3°  to  +3°  will  be  traced.  The  k  value  for 
the  rays  is  the  same  whether  the  angle  is  positive  or  negative* 
When  the  rays  going  up  are  traced*  one  also  gets  the  corresponding 
rays  going  down.  It  is  only  a  matter  of  shifting  the  rays  a 
constant  horizontal  distance  until  the  crossing  of  the  source 
depth  coincides  with  the  source. 

Table  3  shows  the  result  of  the  k  value  program. 


TABLE  3 

V 

1527.13 

S 

1 

S 

0.99984770 

b0 

1527.36261732 

aO 

2 

s 

0.99939083 

bO 

1 528 .06084882 

aO 

3 

s 

0.99862954 

bO 

1529.22574271 

aO 
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Sound  Speed  m/s  Range  km 

1500  1510  1520  1530  0  1  2 


UI  U.>doQ 
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F'j.  6  EXAMPLE  OF  RAY-TRACING  BY  DTC 


We  call  the  start  of  the  program  by  Vj  and  give  initial 
data  ( c Q  =  1527*  13  m/s)  and  press  5*  Then  we  give  the 
data  (  Q0  =  1°)  and  press  S  and  the  DTC  prints  cos  0o(BO) 
and  the  k  value  (AO)*  Thereafter  we  continue  in  the  same 
way  for  other  values  of  0O. 

The  program  is  used  with  many  decimal  places  because  the  cosine 
for  small  angles  is  very  near  to  unity.  But  we  select  only  two 
decimal  places  when  we  use  the  result  as  input  data  in  the  ray 
tracing  program. 

The  computation  of  the  ray  starting  with  a  positive  angle 
of  2°  is  shown  in  Table  4* 

First  we  call  the  start  of  the  program  by  V  and  give  the 
initial  data  (k  =  1528. 06,  cQ  =  1527*13  and  zQ  =  10), 
pressing  5  after  each  initial  data. 

The  DTC  prints  out  tan  0O(A^/  which  is  needed  for  the 
intensity  program . 

Then  we  ask  for  horizontal  distance  to  the  point  where  the 
ray  reaches  the  surface, by  giving  the  sound  speed  at  the 
surface  (1528.37)  and  the  depth  as  input  data.  But  as  the 
sound  speed  at  the  surface  is  bigger  than  the  k  value  of 
the  ray,  the  surface  cannot  be  reached.  The  computer 
therefore  prints  out  the  depth  of  the  turning  point  ( bO ) 
and  the  horizontal  distance  to  the  turning  point  (D^)  . 

The  r  1  (  C  0)  for  the  intensity  program  and  the  tan  0(A^) 
are  also  printed  out. 

Next  to  be  asked  for  is  the  horizontal  distance  to  the  point 
where  the  ray  returns  to  the  source  depth.  As  the  ray  will  now 
cross  this  depth  at  the  same  angle  as  the  initial  angle,  we 
now  also  get  the  data  for  the  ray  starting  with  an  initial 
angle  of  -2° .  The  table  shows  how  the  computation  goes  on 
until  the  ray  is  plotted. 
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TABLE  4 


V 

1506.01 

S 

1520.06 

S 

75 

S 

1527.13 

S 

-15356.1300 

CO 

10 

S 

1412.0645 

DO 

0.03494 

AO 

0.1721 

A  0 

1520.37 

S 

1506.59 

S 

0 

S 

90 

S 

2.5000 

bO 

1429.7094 

CO 

-445.0233 

CO 

1501.7924 

DO 

429.9074 

D  0 

0.1699 

AO 

0.0000 

AO 

1507.93 

S 

1527.13 

S 

160 

S 

10 

S 

2007.0700 

co 

-445.0233 

cO 

'*926.0206 

DO 

059.0140 

DO 

0.1642 

AO 

0.03^4 

AO 

1509.02 

S 

1524.22 

s 

210 

S 

20 

s 

1462.4094 

CO 

-109.9720 

co 

2241. 6250 

DO 

1040.0757 

DO 

o.-’sss 

AO 

0.0722 

AO 

1511.51 

s 

1511.11 

s 

310 

s 

30 

s 

3271.130'’ 

co 

-42.1670 

CO 

2904.6030 

DO 

1139.0527 

DO 

0.1409 

AO 

0.1502 

AO 

1506.10 

s 

50 

s 

•220.6070 

co 

1265.4456 

DO 

0.1715 

AO 

30 


The  plotting  of  the  other  rays  is  done  in  the  same  way,  the 
result  being  shown  in  Fig.  6b.  It  is  easily  seen  from  the 
figure  that  the  ray  that  is  tangent  to  the  surface  is  the 
limit  of  the  shadow  zone. 

As  we  have  no  crossing  of  rays  in  this  case,  the  travel-time 
computation  is  of  little  interest.  The  intensity  calculation, 
however,  is  of  great  interest,  because  Fig.  6b  clearly  shows 
that  the  intensity  in  the  domain  covered  by  rays  with  a 
positive  initial  angle  is  less  than  that  for  negative  angles. 

Table  5  shows  intensity  calculation  for  two  rays.:  those  with 
initial  angles  of  -2°  and  +2°  .  For  the  ray  with  -2° 
initial  angle,  the  intensity  is  asked  for  at  distances  of 
1266  m  and  1926  m.  By  taking  20  log  of  the  answer  (AO*) 
one  finds  the  geometric  spreading  losses  to  be  66.6  dB  and 
68.1  dB  respectively.  By  taking  20  log  of  the  next  answer, 

(AO)  =  (Vlo/lp)  /X,  one  obtains  a  comparison  with  the  spherical 
losses.  In  this  case  it  turns  out  that  the  geometrical 
spreading  losses  are  respectively  1.6  dB  and  2.4  dB  higher 
than  the  spherical  spreading  losses.  For  the  ray  that  starts 
with  initial  angle  +2°  the  intensity  is  computed  at  a  distance 
of  1066  m,  and  the  intensity  here  is  2.2  dB  higher  than 
spherical  spreading. 

The  same  data  have  been  run  with  SACLANTCEN’s  Ray  Tracing 
Program  (C.P.  67)  on  SACLANTCEN T s  Elliott  503  computer. 

The  result  for  one  initial  angle  is  shown  in  Table  6. 

Table  7  compares  the  results  of  the  two  ray— tracing  programs 
for  the  -2°  initial  angle.  Under  the  travel-time  one 
column  shows  the  correction  for  the  difference  in  path  length. 

3 . 2  Use  of  a  pre-computation  to  prepare  ray-tracing 

data  for  a  large  computer 


The  above  ray-tracing  example  was  to  illustrate  how  to  use 
the  DTC  for  normal  ray  tracing.  Nowadays  at  most  of  the  places 


TABLE  5 


-2°  initial  angle 


V 

1528.06  S 

1527.13  S 

0.03494  S 

1528.37  S 

-445.8  S 

1527.13  S 

-445.8  S 

1524.22  S 

-190  S 

1511.11  S 

-42.2  S 

1506.1  S 

-220.7  S 

Z 

1265  S 

2149.68090  A  0 

1.69935  A  Q 

1506.01  S 

-15356.2  S 

1506.59  S 

1429.7  S 

1507.93  S 

2887.9  S 

Z 

1926  S 

2544.63568  A  <) 

1.32120  A  0 


+2°  initial  angle 


V 

1528.06  S 
1527.13  S 
0.03494  S 

1524.22  S 
-190  S 

1511.11  S 
-42.2  S 

1506.1  S 
-220.7  S 

1506.01  S 
-15356  S 

1506.59  S 
1429.7  S 

1507.93  S 
2887.9  S 

Z 

1066  S 
830.18232  A  0 
0.77878  A  0 
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TABLE  6 


Raytr, A.S.TM 

Initial  angle  -1.999 

CODE  0=D  1-U  2 

=S  3=*B 

4-UD 

5=DU 

CODE 

distance 

depth 

travel- 

angle 

path- 

intensity 

time 

length 

m 

m 

sec 

degrees 

m 

dB 

4 

429.9 

2.5 

.281422 

0.00 

429.9 

-52.67 

0 

859.7 

10.0 

.562846 

2.00 

859.9 

-58.69 

0 

1048.6 

20.0 

.686802 

4.06 

1049.0 

-62.13 

0 

1139.1 

30.0 

.7468'' 1 

8.54 

1140.1 

-65.65 

0 

1263.5 

50.0 

.830327 

9.73 

1266.1 

-66.62 

0 

1409,2 

75.0 

.928492 

9.75 

1414.0 

-67.07 

0 

1497.2 

90.0 

.987712 

9.62 

1503.2 

-67.25 

0 

1917.1 

160.0 

1.27018 

9.31 

1928.9 

-68.08 

0 

2229.8 

210.0 

1.48002 

8.86 

2245.6 

-68.43 

0 

3500.0 

398.8 

2.32967 

8.05 

3529.7 

-69.50 

where  one  does  ray  tracing,  one  normally  has  access  to  a  computer. 
But  even  then  the  DTC  and  the  programs  described  here  can  be  a 
help.  For  a  certain  sound  speed  gradient,  one  can  make  a 
short  and  fast  ray  tracing  by  means  of  the  DTC.  Afterwards, 
it  is  much  easier  to  ask  for  the  correct  things  from  the  large 
computer . 

The  example  given  in  Fig.  7  will  clarify  this.  The  sound  speed 
profile  is  shown  at  the  left.  On  the  right  are  some  rays 
plotted  by  means  of  the  DTC  ray-tracing  program.  The  rays 
selected  are  those  that  have  turning  points  at  depths  where 
the  sound  speed  gradient  has  discontinuity  points.  This 
simple  ray  tracing  can  give  a  good  idea  about  shadow  zones, 
focusing,  multipath  area,  or  caustics,  thereby  allowing  more 
exact  data  to  be  specified  when  using  a  large  computer. 
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COMPARISON  BETWEEN  THE  RESULTS  OF  RAY-TRACING  PROGRAMS 
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PREPARATION  op  RAY  TRACING  DATA  FOR  USE  With  a  large  computer 
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FIG.  8  TRAVEL-TIME  ANOMALY  AS  A  FUNCTION  OF  DISTANCE  AT  600  m  DEPTH  FOR  DIFFERENT  iNiTiAL  ANGLES 


A  comparison  was  made  between  the  results  of  these  different 
DTC  programs  and  those  given  by  the  programs  of  a  large 
computer,  using  the  same  sound  speed  profile  [Fig.  6].  In 
spite  of  the  comparatively  extreme  conditions  shown  by  this 
profile  —  thick  layers,  rather  large  sound  speed  differences, 
and  initial  angles  of  up  to  10°  —  the  comparison  revealed 
that  the  DTC  program  makes  an  error  of  less  than  1%  in 
horizontal  range  and  of  less  than  in  travel  time 

(corrected  for  difference  in  horizontal  range).  By 
dividing  into  thinner  layers, the  error  in  corrected  travel 
time  can  be  reduced  to  less  than  0.3^°*  The  error  in 
intensity  program  is  normally  less  than  0.1  dB ,  but  in 
cases  with  intensity  anomalies  of  the  order  of  20  dB, 
the  error  is  less  than  1  dB * 


3 . 3  Example  of  Travel -time  Program 

From  Fig.  7  we  see  that  rays  with  initial  angles  of  about 
0°  cross  rays  with  initial  angles  of  about  -10°.  To 
analyze  this  in  more  detail,  some  rays  have  been  computed 
by  means  of  the  ray-tracing  program  and  travel— time  program. 
(The  layer  between  100  m  and  3 5C  m  was  divided  into  four 
layers  during  the  computation  to  get  good  accuracy.) 

In  Fig.  8  the  travel-time  anomaly  is  plotted  as  a  function 
of  distance  at  600  m  depth  for  different  initial  angles.  We 
notice  that  when  the  initial  angle  increases  from  0°,  the 
horizontal  distance  out  to  a  depth  crossing  of  600  m  first 
decreases  up  to  about  1°  initial  angle  and  afterwards 
increases.  For  the  other  family  of  rays,  we  know  that  for 
initial  angles  steeper  than  -10.05°  surface  reflection 
occurs  and  the  horizontal  distance  decreases  sharply. 
Therefore,  it  is  only  at  about  12040  m  that  intersection 
of  the  two  ray-families  can  occur  at  600  m  depth  (neglecting 
short  distance  with  surface-reflected  r»ays)  . 

We  see  from  the  figure  that  the  ray  with  about  -10°  initial 
angle  arrives  about  39  ms  before  the  one  with  about  1° 
angle . 
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initial 


CONCLUSION 


By  means  of  the  programs  presented  here  and  a  DTC,  it  is 
possible  to  compute  the  sound  field  in  the  same  manner  as  with 
a  ray-tracing  program  for  a  large  computer.  Of  course  the 
time  necessary  for  the  computation  is  quite  different. 

For  a  large  program  the  Olivetti  101  Desk  Top  Computer  uses 
about  15  seconds  for  each  sequence.  It  should  therefore 
take  about  one  hour  to  convert  a  set  of  temperature,  salinity, 
and  depth  profiles  into  a  sound  speed  profile,  make  ray-tracing 
for  some  rays,  and  make  intensity  calculation  at  some  points 
(assuming  about  ten  layers). 

Where  a  computer  is  not  available,  it  is  thus  possible  by 
means  of  a  DTC  to  make  some  ray-tracing  within  a  reasonable 
time.  When  a  computer  is  available  but  is  heavily  booked, 
the  use  of  these  programs  on  a  DTC  can  be  useful  to  give 
some  results  while  waiting  to  obtain  the  use  of  the  large 
computer . 
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TABLE  9 


LEROY  2m  SOUND  SPEED  FORMULA 


REGISTER  1 

REGISTER  2 

REGISTER  F 

REGISTER  E 

REGISTER  D 

1 

A  V 

25 

RO 

49 

C  t 

73 

B  ♦ 

97 

2 

A/t 

26 

R/{ 

50 

C/  : 

74 

0/  A 

98 

3 

R/  S 

27 

D*' 

51 

B  + 

75 

A/,  t 

99 

4 

D  ^ 

28 

X 

52 

B* 

76 

R  - 

100 

5 

C/t 

29 

B  + 

53 

B/  A 

77 

0/  s 

101 

6 

S 

30 

. 

54 

A/  t 

78 

- 
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7 

X 

31 

B/  x 

55 

R  0 

79 

c/t 
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8 

C/  : 

32 

Ct 

56 

0/  A 

80 

c/  A 
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9 

B/$ 

33 

c  - 

57 

81 

A/  t 

105 

10 

B/V 

34 

Bt 

58 

C  X 

82 

0/  t 
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11 

A/t 

35 

5 

59 

C/  : 

83 

- 
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12 

F  x 

36 

C  r 

60 

A  x 

84 

B/  x 

108 

13 

E/  + 

37 

c  * 

61 

A  + 

85 

C/  x 
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14 

X 

38 

A/t 

62 

B  + 

86 
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R/a 

63 
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40 

D  x 

64 

S 

88 

R  t 

112 

17 

D/  + 

41 

: 

65 

X 

89 

0/  s 

113 

18 

M 
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90 

X 
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19 

S 

43 

Bt 

67 

R/  S 

91 
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20 

i 

44 

c  x 

68 

D  * 

92 

A  6 
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21 

c/  : 

45 

C/  : 

69 

: 

93 

/O 

117 

22 

46 

: 

70 

C  x 

94 

V 
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23 

c  * 

47 

; 

71 

c/  : 

95 

119 

24 

A/t 

48 

A  x 

72 

B  + 

96 

120 

CONSTANTS 


CONSTANTS 


1401.14 


0/ 


49.1 
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TABLE  10 


SOUND  SPEED  -  DEEP  OCEAN  TERM  AND  LOW  SALINITY  TERM 


REGISTER  1 

REGISTER  2 

REGISTER  F 

REGISTER  E 

REGISTER  D 

1 

A  V 

25 

8  T 

49 

73 

97 

2 

S 

26 

50 

74 

98 

3 

l 

27 

0/  4- 

51 

75 

99 

4 

0/  : 

28 

s 

52 

76 

100 

5 

A  x 

29 

+ 

53 

77 

101 

6 

8  T 

30 
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54 

78 
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7 

31 

C/  x 

55 

79 

103 

8 

F  + 

32 

F/  x 

56 

80 
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9 

0/ 1 

33 

D/  : 

57 

81 
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10 

B/  J- 

34 

58 

82 

106 

11 
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35 
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59 

83 
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12 

B/  x 

36 

E/  - 

60 

84 
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13 
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37 

8  x 

61 

85 
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14 
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38 
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62 

86 
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15 
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63 

87 
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64 

88 
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21 
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93 
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70 

94 

118 

23 

D  - 

47 

71 

95 
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24 

A  x 

48 

72 

96 

120 
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0/  t 

3 

B  t 

35 

D  t 

1.5 

F/  1 

30 

E/  t 

4 

F  t 

43 


TABLE  11 


K -VALUE  PROGRAM 


REGISTER  1 

REGISTER  2 

REGISTER  F 

REGISTER  E 

REGISTER  D 

1 

A  V 

25 

49 

73 

97 

2 

S 

26 

50 

74 

98 

3 

C  T 

27 

51 

75 

99 

4 

A  W 

28 

52 

76 

100 

5 

S 

29 

53 

77 

101 

6 

i 

30 

54 

78 

102 

7 

F  : 

31 

55 

79 

103 

6 

A  x 

32 

56 

60 

104 

9 

F/  : 

33 

57 

81 

105 

10 

A  X 

34 

58 

82 

106 

11 

B  t 

35 

59 

83 

107 

12 

E  : 

36 

60 

84 

108 

13 

B  + 

37 

61 

85 

109 

14 

E/  + 

38 

62 

86 

110 

15 

B  | 

39 

63 

87 

111 

16 

B  0 

40 

64 

88 

112 

17 

c  * 

41 

65 

89 

113 

18 

B  : 

42 

66 

90 

114 

19 

A  0 

43 

67 

91 
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20 

to 

44 

68 

92 
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21 

W 

45 

69 

93 
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22 

V 

46 

70 

94 
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23 

47 

71 

95 
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24 

48 

72 

96 
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CONSTANTS 

1 

E/  t 
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E  t 

6 

E  t 

I 

-2 
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TABLE  12 


RAY-TRACING  PROGRAM 


REGISTER  1 

REGISTER  2 

REGISTER  F 

REGISTER  E 

REGISTER  D 

1 

A  V 

25 

C/i 

49 

c/  + 

73 

M 

97 

2 

S 

26 

c  t 

50 

B  - 

74 

C  0 

98 

3 

B  T 

27 

A/  V 

51 

C  1 

75 

D  + 

99 

4 

S 

28 

s 

52 

B  ± 

76 

0  * 

100 

5 

C/  T 

29 

l 

53 

C/* 

77 

D  0 

101 

6 

s 

30 

B/  ♦ 

54 

B/0 

78 

A  Z 

102 

7 

B/  t 

31 

B/  - 

55 

A  Y 

79 

c/  i 

103 

8 

0  * 

32 

□/♦ 

56 

B  * 

80 

B  : 

104 

9 

z 

33 

C/  * 

57 

A  + 

81 

□/$ 

105 

10 

A  W 

34 

B  - 

58 

A 

82 

0/  * 

1C5 

11 

s 

35 

/  w 

59 

D/  x 

83 

A  x 

107 

12 

36 

Y 

60 

C  : 

84 

A  : 

108 

13 

37 

A/  W 

61 

85 

- 

109 

14 

C/  - 

38 

A  - 

62 

A  - 

86 

A  J~ 

110 

15 

c  ♦ 

39 

- 

63 

- 

87 

D  /  : 

111 

IS 

C  * 

40 

D/  x 

64 

C/  - 

88 

aO 
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17 

A  ♦ 

41 

C  : 

65 

B  + 

89 

/o 

113 

18 

/  v 

42 

B/$ 

66 

a  r 

90 

W 

114 

19 

C/  * 

43 

B/  + 

67 

91 

115 

20 

A/  t 

44 

B/t 

68 

B  * 

92 

116 

21 

F  1 

45 

□/$ 

69 

C/  - 

93 

117 

22 

F  S 

46 

D/  - 

70 

a  r 

94 

118 

23 

E/  S 

47 

D/t 

71 

0/  - 

95 

119 

24 

48 

c  x 

72 

C  X 

96 
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CONSTANTS 

t 

t 

t  t 

f  t 

! 

t 
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TABLE  13 


INTENSITY  PROGRAM 


REGISTER  1 

REGISTER  2 

REGISTER  F 

REGISTER  E 

REGISTER  D 

1 

A  V 

25 

At 

49 

t 

73 

8/  x 

97 

2 

C  * 

26 

/  V 

50 

D  ± 

74 

A  0 

90 

3 

S 

27 

D/  A 

51 

: 

75 

8/  : 

99 

4 

B  T 

20 

A/t 

52 

76 

100 

5 

S 

29 

R  4 

53 

E/  - 

77 

A  0 

101 

6 

□/  r 

30 

R  S 

54 

s 

78 

/O 

102 

7 

B  4. 

31 

E/  S 

55 

X 

79 

w 

103 

0 

0/  : 

32 

56 

C/  X 

80 

104 

9 

s 

33 

D/t 

57 

c  + 

01 

105 

10 

X 

34 

A/  V 

50 

ct 

82 

106 

11 

c/t 

35 

B  ± 

59 

/o 

03 

107 

12 

B  4 

36 

0/  - 

60 

W 

84 

100 

13 

□ /  - 

37 

/  w 

61 

A  2 

05 

109 

14 

A  V 

30 

D  * 

62 

S 

06 

110 

15 

> 

39 

A  - 

63 

r 

07 

1  Vi 

16 

D  4 

40 

- 

64 

B/t 

00 

112 

17 

41 

65 

B  4 

89 

113 

18 

E/t 

42 

s 

66 

A/  t 

90 

114 

19 

/o 

43 

/o 

67 

D/  T 

91 

115 

20 

A  W 

44 

w 

68 

: 

92 

116 

21 

S 

45 

A/  W 

69 

A^ 

93 

117 

22 

4 

46 

B  * 

70 

E/  : 

94 

110 

23 

D/i 

47 

0/  - 

71 

C  x 

95 

119 

24 

D/  - 

40 

a  r 

72 

A  J~ 

96 

120 

CONSTANTS 


CONSTANTS 


46 
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TABLE  14 


TRAVEL-TIME  PROGRAM 


REGISTER  1 

REGISTER  2 

REGISTER  F 

REGISTER  E 

REGISTER  D 

1 

A  V 

25 

S 

49 

w 

73 

97 

2 

S 

26 

- 

50 

74 

96 

3 

8  T 

27 

0/  t 

51 

75 

99 

4 

s 

28 

M 

52 

76 

100 

5 

C/  T 

29 

M 

53 

77 

101 

6 

C  t 

30 

B  * 

54 

78 

102 

7 

D  * 

31 

8/  - 

55 

79 

103 

8 

0/  * 

32 

E  x 

56 

80 

104 

9 

/O 

33 

E/  x 

57 

81 

105 

10 

A  W 

34 

8  : 

58 

82 

106 

11 

S 

35 

* 

59 

83 

107 

12 

l 

36 

D  + 

60 

84 

108 

13 

8  - 

37 

DJ 

61 

85 

109 

14 

/  V 

38 

D/X 

62 

66 

110 

15 

B  + 

39 

B  : 

63 

87 

111 

16 

V 

40 

D  + 

64 

68 

112 

17 

A/  V 

41 

A  0 

65 

69 

113 

18 

B  ^ 

42 

B/$ 

66 

90 

114 

19 

A  Y 

43 

D/  x 

67 

91 

115 

20 

44 

C/  : 

68 

92 

116 

21 

C  + 

45 

B/J 

69 

93 

117 

22 

E/  : 

46 

B/  - 

70 

94 

ns 

23 

B/$ 

47 

A  0 

71 

95 

119 

24 

□/  4 

48 

/o 

72 

96 

120 

CONSTANTS 

CONSTANTS 

2 

E/  f 

t 

t 

f  t 

t 

t 

47 


APPENDIX  A 


APPLICATION  OF  PROGRAMS 


TO  AN  OLIVETTI  101  DESK-TOP  COMPUTER 


APPENDIX  A 


The  appendix  contains  first  a  brief  description  of  an  Olivetti  101 
Desk  Top  Computer,  enabling  the  reader  to  follow  the  programming. 

Then  follows  an  easy  program  (pressure-depth  program)  in  which 
every  instruction  in  the  programming  is  explained. 

The  last  part  contains  some  remarks  about  the  application  of 
the  formulae  to  an  Olivetti  101  together  with  the  complete 
programs  in  Tables.  Instructions  for  the  use  of  an  Olivetti  101 
are  also  given. 


A. 1  THE  OLIVETTI  101  ELECTRONIC  DESK-TOP  COMPUTER 

The  desk-top  computer  works  with  programs  recorded  on  magnetic 
cards.  As  shown  in  Fig.  A.l  (taken  from  Ref.  8)  it  consists  of 
the  following  basic  units: 

a.  A  Memory  to  hold  numerical  data  and  program 
instructions . 

b.  A  Keyboard  to  put  data  or  instructions  into 
the  computer. 

c.  A  Printing  Unit. 

d.  A  Read  or  Recording  Unit  for  the  magnetic 
ca  rds . 

e.  A  Control  and  Arithmetic  Unit  to  receive  instructions 
from  the  memory,  interpret  their  meeting,  and  carry 
them  out  through  the  arithmetic  unit. 
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FiG.  A  I  BLOCK  DIAGRAM  OF  THE  OLIVETTI  101  ELECTRONIC  DESK-TOP  COMPUTER 

The  memory  consists  of  ten  registers  —  two  program  instruction 
registers,  three  operational  registers,  and  five  storage 
registers.  Each  of  the  storage  registers  can  be  split  into 
two  parts.  The  two  program  instruction  registers  can  store 
48  instructions.  If  more  than  this  number  are  given,  he 
computer  starts  to  fill  one  of  the  five  storage  registers 
with  instructions  and  it  can  then  no  longer  be  used  for 
storage.  Up  to  three  of  the  storage  registers  can  be 
filled  with  instructions,  thus  enabling  the  computer  to  take 
120  instructions. 

The  computer  is  capable  of  normal  arithmetical  operations  and 
in  addition  can: 

a.  Give  absolute  values. 

b.  Extract  square  root  of  absolute  values 

c.  Generate  constants  by  means  of  pseudo¬ 
program  instruct! ons . 
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d.  Make  conditional  and  unconditional  jumps  in 
the  instruction  sequence. 

The  content  of  the  two  program  instruction  registers  and  three 
of  the  storage  registers  can  be  recorded  on  a  magnetic  card, 
and  afterwards  read  into  the  computer  by  the  same  card. 

This  means  that  the  "program  card"  can  contain  only  program 
instructions,  or  can  contain  program  instructions  and  constants. 

A. 2  PROGRAM  FOR  THE  USE  OF  A  PRESSURE  FORMULA 

A  simple  formula  [Ref.  9]*  which  gives  the  absolute  pressure 
in  the  ocean  as  a  function  of  depth  and  latitude,  has  been 
proposed  by  C.C.  Leroy.  This  formula  has  sufficient  accuracy 
for  use  in  underwater  acoustics. 

The  formula  is: 

P(z)  =  1 .04  +  o.  102506  (  1  +  0.00528  sin2  cp)  Z  +  2  .  524  X  10_7  z2  , 

where  P  is  in  kg/cm2  ,  z  is  in  metres,  and  cp  is  the 
latitude  in  degrees. 

To  compute  this  formula  on  a  DTC,  it  is  first  necessary  to  make 
a  program  for  sin  cp.  This  is  done  by  computing  sin  cp  after 
the  mathematical  series 


X  is  the  corresponding  value  in  radians. 

In  the  program,  four  terms  of  the  series  are  used.  The  pressure 
formula  may  be  written  as  follows: 

P(z)  =  1.04+  10.2  506(  10  +  0.0  528  sin2  cp)  z/l03+  0.2  524  z2/l06. 
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The  accuracy  of  the  sin  ep  part  of  the  program  is  best  for  small 
angles,  but  even  for  90°  the  result  is  0.99992. 

The  program  is  shown  in  Table  8  of  the  main  text  and  will  be 
explained  in  detail. 

The  program  consists  of  57  instructions,  which  means  that  two 
program  registers  and  half  the  F  register  are  occupi'cLby  program 
instructions.  This  leaves  space  for  five  constants.  These 
constants  are  put  into  the  memory  after  the  program  instructions 
are  made,  and  afterwards  both  the  program  instructions  and 
the  constants  are  recorded  on  the  magnetic  card. 

For  convenience,  the  constants  used  are  each  allocated  a 
letter : 


k  =  1000,  m  =  57.29578,  n  =  10.2506,  p  =  0.0528,  q  =  0.2524. 

It  is  important  to  note  that  the  arithmetical  operations  take 
place  with  the  A  and  M  register,  so  the  content  of  a  register 
is  first  transferred  into  the  M  register  before  an  arithmetical 
operation  takes  piace  together  with  the  A  register. 

An  explanation  of  the  various  instructions  follows: 

1.  Label  at  beginning  of  the  program.  N.B.  the  constants 
are  automatic  in  their  registers,  when  the  program 
starts . 

2.  Stop  to  receive  input  from  keyboard. 

3.  Put  input  data  (z  =  depth)  into  A  register. 

4.  Divide  content  of  A  register  by  content  of  d  register. 

5.  Exchange  content  of  A  register  with  content  of  C 

register  (which  means  here:  put  content  of  A  register 

into  C  register) . 
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6.  Stop  to  receive  input  from  keyboard. 

7.  Put  input  data  (latitude  in  degrees)  into  A  register. 

8.  Take  absolute  value  of  content  of  A  register. 

9.  Conditional  jump.  If  content  of  A  register  is  negative 
or  zero,  the  computer  goes  on  in  the  sequence;  if  it  is 
positive  the  computer  jumps  to  A/V.  As  the  absolute 
value  of  the  content  in  the  A  register  is  taken,  it 
goes  on  in  the  sequence  only  if  the  input  latitude 

is  zero. 

10.  Means  go  to  label  AW. 

1  1  .  See  9  • 

12.  Divide  content  of  A  by  content  of  D;  convert  degrees 
to  radians. 

13 •  Multiply  content  of  A  by  itself  and  therefore  also 

transfer  the  content  of  A  (before  the  multiplication) 
into  M. 

14*  Multiply  A  and  M,  and  now  we  have  X  in  the  A. 

15*  Transfer  content  of  M  into  B. 

16.  and  17*  Generate  the  figure  -6  in  M, 

18 .  Divide  A  by  M,  and  we  now  have  -X  /6  =  -X  /3  l  in  A. 

19*  Multiply  content  of  A  register  by  itself',  also 

transferring  content  of  A  (before  multiplication)  into  M, 

20.  Transfer  content  of  M  (-X3/6)  into  b. 

Divide  content  of  A(X^/36)  by  content  of  B. 
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21. 


22,  23 ,  24  and  25*  Generate  the  figure  3*33  in  M. 

26.  Divide  content  of  A  by  content  of  M,  making  content 

of  A  equal  to  X^/51  • 

27.  Exchange  content  of  A  with  content  of  c. 

28.  Transfer  content  of  c  into  A  (also  leaving  content 
of  c  unchanged) . 

29.  Multiply  content  of  A  by  content  of  b. 

30.  Divide  content  of  A  by  content  of  B. 

31.  and  32.  Generate  figure  7  in  M. 

33*  Divide  content  of  A  by  content  of  M.  Now  the 

7 

content  of  A  is  -X/7*  • 

34*  Add  content  of  c  (X^/50  to  content  of  A. 

35*  Add  content  of  b  (-X3/3*)  to  content  of  A, 

36.  Add  content  of  B(X)  to  content  of  A.  Now  the 
content  of  A  equal  X-X3/3J  +  X 5/ 5 1  -  X7/7!  . 

37*  Label  which  the  computer  jumps  to  in  case  the  input 
data  equal  0°  (see  10) . 

38.  Multiply  content  of  A  by  itself,  making  sin2  cp  in  A 

39*  Multiply  contents  of  E  with  contents  of  A,  making 
0.0528  sin2  cp  in  A. 

40,  4  1>  42.  Generate  figure  10  in  M. 

43*  Add  content  of  M  to  content  of  A  making  10  +  0.0528 
in  A . 


sin2  cp 
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44*  Multiply  content  of  C  by  content  of  A,  making 
(  10  +  0.0528  sin2  cp)  Z/1000  in  A. 

45.  Multiply  content  of  e  by  content  of  A,  making 
10.2506(  10  +  0.0528  sin2  cp)  Z/1000  in  A. 

46.  Exchange  content  of  C  with  contents  of  A,  which 
means  that  the  result  in  A  is  put  into  C,  and  the 
content  of  C  (Z/lOOO)  is  put  into  A. 

47*  Multiply  the  content  of  A  with  itself  making 
Z2 /  10^  in  A . 

48.  Multiply  content  of  f  with  content  of  A,  making 

0.2524  Z2/l06. 

49 .  Add  content  of  A  to  content  of  C,  making 

10.2506  (  10  +  0.0528  sin2  cp)  Z/1000  +  0. 2524  Z2/10°  in  A 

50.  51 5  52,  53*  Generate  the  figure  1 . 04  in  M. 

54*  Add  the  content  of  M  to  the  content  of  A,  completing 
the  formula  in  A. 

55 .  Print  out  result  of  A. 

56.  Advance  the  paper  without  printing. 

57*  Instruction  to  go  to  label  AV  (at  beginning  of  the 
program) . 


Constants  to  be 

1000 

57.29578 
10. 2506 
0.0528 
0.2524 


put  into  register: 

into  register  d 
into  register  D 
into  register  e 
into  register  E 
into  register  f 
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Instruction  for  use  of  the  program 


a.  Read  the  program  into  the  DTC  from  the  magnetic 

card  (the  knob  Reg.  Progr.  must  be  at  its  outermost  position), 

b.  Put  the  decimal  selector  on  5 »  This  is  because 
of  the  accuracy  during  the  computation,  and  not  because  we 
want  there  to  be  this  number  of  decimal  places  in  the  result. 

c.  Call  program  by  pressing  V. 

d.  Print  z  (depth  in  metres)  on  keyboard  and  press  S. 

e.  Print  cp  ( latitude  in  degrees)  on  keyboard  and 

press  S. 

f.  Read  result  (pressure  in  kg/cm2 )  printed  out. 

g.  Give  new  data. 

A. 3  SOUND  SPEED  PROGRAMS 

The  programs  for  the  main  formula  contain  a  lot  of  instructions, 
so  there  are  only  two  half  registers  for  constants:  1401.14 

in  d  and  49  1  in  D,  The  program  is  shown  in  Table  9  and 

the  flowchart  in  Fig.  2  of  the  main  text. 

The  program  for  Vc  and  Vd  has  less  than  48  instructions, 
so  there  is  space  for  six  constants  in  the  registers. 

These  are 


-1000 


d 


35 

30 

3 


D 


E 


e 


1.5 

4 


F 


f 


The  program  is  shown  in  Table  10  of  the  main  text. 
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Instruction  for  Use 


Main  Sound  Speed  Program 

Read  the  program  into  the  DTC.  Put  decimal  selector  on  4 
and  call  program  by  V. 

Input  data: 

t  :  temperature  in  degrees  centigrade, 

s  :  salinity  in  %0  . 
z  :  depth  in  metres, 
cp  :  latitude  in  degrees. 

Read  sound  speed  in  m/s 
Give  new  data. 

Low  salinity  or  deep  ocean  terms 

Put  decimal  selector  on  3  «> 

If  Vc  (deep  ocean)  is  wanted,  press  V  and  give  depth  in  m 
as  input  data,  read  Vc  in  m/s.  Vc  must  then  be  added  to  the 
sound  speed  calculated  by  means  of  the  program  for  the  basic 
formula.  If  Vd  (low  salinity  term)  is  wanted,  press  V\ , 

Give  temperature  (°C),  salinity  (%>)  ,  depth  (m)  as  input, 
read  Vc  (m/s)  as  output.  Vc  must  then  be  added  to  the 
basic  formula. 


A. 4  PROGRAM  TO  COMPUTE  k-VALUE 

The  program  is  shown  in  Table  11  of  the  main  text. 

Read  the  program  into  the  DTC,  put  decimal  selector  on  8 
and  call  program  by  V. 


56 


Initial  data: 


cQ:  sound  speed  at  source  in  m/s. 


Data  : 

0O:  initial  angle  in  degrees. 

Result : 

co 

cos  0.  and  k  =  - —  . 

0  cos  0O 


New  data: 

For  new  initial  data,  press  V. 


A. 5  RAY  TRACING  PROGRAM 

The  flowchart  for  the  program  is  shown  in  Fig.  3  and  the 
program  in  Table  12  of  the  main  text. 

Read  the  program  into  the  DTC. 

Decimal  selector:  during  initial  data  computation  on  5> 
otherwise  on  4* 

Call  program  by  V. 

Initial  data  : 

c 

k  :  ray  constant  =  - r —  , 

cos  0^ 

cQ  :  sound  speed  at  source  m/s y 

zQ  :  depth  of  source  in  m  . 

Read  tan  0 Q ,  (  @Q  initial  angle).  This  result  is  needed  for 

the  intensity  program. 
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Data  : 


c^:  sound  speed  in  m/s  in  layer  n  (n  =  1,  2 . ), 

z^:  depth  of  layer  n. 

Read 

r^  needed  for  intensity  program, 

X  horizontal  distance  in  m, 

n  7 

tan  0  ,  where  0  is  angle  out  of  layer  n. 
n7  n  J 

Data  for  next  layer, 

If  a  sound  speed  is  given  that  is  bigger  than  the  k  value 
of  the  ray ,  a  turning  point  is  passed*  The  DTC  then  prints 
out  the  depth  of  the  turning  point  and  modifies  the  input 
data  * 

After  a  turning  point  or  a  reflection,  one  has  to  trace 
the  ray  back,  giving  once  more  the  depths  and  their 
corresponding  sound  speeds  as  data. 

If  the  operator  does  not  notice  that  a  turning  point  is 
passed  and  asks  for  a  new  point  with  a  sound  speed  higher 
than  the  k  value  of  the  ray,  then  the  DTC  stops  the 
computation  and  the  red  light  comes  on.  The  main  loop  of 
the  program  can  be  called  by  W,  and  the  data  stored  in  the 
computer  are  still  correct. 

When  a  ray  is  traced,  press  V  and  give  initial  data  for 
next  ray. 

A. 6  INTENSITY  PROGRAM 

The  flowchart  of  the  program  is  given  in  Fig.  4  and  the 
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program  in  Table  13  of  the  main  text.  Constants  to  be 
put  in:  1  into  D. 

The  program  is  based  on  the  ray  tracing  program  being  used 
first,  and  the  values  for  k,  tan  0O  ,  and  r^  are  therefore 
known . 

Load  the  program  by  means  of  the  magnetic  card,  and  put 
decimal  selector  on  5* 

Call  program  by  means  of  V. 

Give  initial  data 


k  constant  for  the  ray, 

cQ  sound  speed  at  the  source, 

tan  0O,  where  0Q  is  initial  angle. 

Thereafter  give 


c^  sound  speed  in  layer  n  (n  =  1,  2 . ), 

r^T  modified  curvature  radius  printed  out  from  ray 

tracing  program  (the  sign  is  important)  until 

layer  p,  where  one  wants  the  intensity,  is 

reached,  and  until  c  and  r  T  are  given  as 

p  p 

input  data , 


Then  press  Z  and  give  ( horizontal  distance  to  the 

point  where  the  ray  leaves  the  layer  p) . 


The  outputs  are 


and 


If  more  intensity  points  are  wanted,  continue  by  giving 
c  .  „  and  r  T.  .  as  input  data,  and  so  on. 

p+1  p+1 

For  limitation  of  the  program,  see  Sect.  2.4  of  main  text. 
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A, 7  TRAVEL  TIME  PROGRAM 


The  flowchart  for  the  program  is  shown  in  Fig*  5  and  the 
program  in  Table  14  of  the  mam  text.  Constants  to  be  put 
into  registers  2  into  e» 

Read  the  program  into  the  DTC„ 


Decimal  selector  on  4,  call  program  by  V. 
Initial  datas 


k  :  constant  for  the  ray, 

cQ  s  sound  speed  at  source  (m/s)* 

Data  s 

c^  s  sound  speed  in  m/s  in  layer  n,  (n  =  1,  2,  „  •  .  »  •  ) , 

:  corresponding  horizontal  distance  (m)  given  by 

ray  tracing  program,. 

Result  2 

t  s  travel  time  in  seconds, 

X 

t  - —  s  travel  time  anomaly  • 
co 


For  new  initial  data,  press  V. 

If  c^  is  given  bigger  than  k,  this  means  that  a  turning 
point  is  reached  and  the  computer  substitutes  c^  by  k„ 
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